Indonesian Districst Getis Filtering

Spatial filtering by Getis

Data, data description & plots

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(spdep) # install.packages("spdep")

Loading required package: sp
Loading required package: spData
To access larger datasets in this package, install the
spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Loading required package: sf
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.4
v tibble  3.0.1     v dplyr   0.8.5
v tidyr   1.0.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(summarytools)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'pryr':
  method      from
  print.bytes Rcpp
For best results, restart R session and update pander using devtools:: or remotes::install_github('rapporter/pander')

Attaching package: 㤼㸱summarytools㤼㸲

The following object is masked from 㤼㸱package:tibble㤼㸲:

    view

library(xtable)

Attaching package: 㤼㸱xtable㤼㸲

The following objects are masked from 㤼㸱package:summarytools㤼㸲:

    label, label<-

library(knitr)
library(ExPanDaR)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(plotly)

Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

library(REAT)

Attaching package: 㤼㸱REAT㤼㸲

The following object is masked from 㤼㸱package:readr㤼㸲:

    spec

library(hdrcde)

This is hdrcde 3.3

library(pdfCluster)

pdfCluster 1.0-3

Attaching package: 㤼㸱pdfCluster㤼㸲

The following object is masked from 㤼㸱package:plotly㤼㸲:

    groups

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    groups

library(readxl)
library(statip)

Attaching package: 㤼㸱statip㤼㸲

The following object is masked from 㤼㸱package:pdfCluster㤼㸲:

    plot

The following object is masked from 㤼㸱package:REAT㤼㸲:

    cv

The following object is masked from 㤼㸱package:sp㤼㸲:

    plot

options(prompt="R> ", digits=3, scipen=999)

library(sf)

cpi <- st_read("Shape file CPI/82 cpi foodstuffs point.shp")

Reading layer `82 cpi foodstuffs point' from data source `D:\Mygithub\tutoring_Felipe\Harry\Shape file CPI\82 cpi foodstuffs point.shp' using driver `ESRI Shapefile'
Simple feature collection with 82 features and 154 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 95.3 ymin: -10.2 xmax: 141 ymax: 5.56
geographic CRS: WGS 84

Warning message:
In readChar(file, size, TRUE) : truncating string with embedded nuls

#cpi <- st_read("Shape file CPI/82 cpi procfood point.shp")
#cpi <- st_read("Shape file CPI/82 cpi housing point.shp")
#cpi <- st_read("Shape file CPI/82 cpi health point.shp")
#cpi <- st_read("Shape file CPI/82 cpi clothing point.shp")
#cpi <- st_read("Shape file CPI/82 cpi education point.shp")
#cpi <- st_read("Shape file CPI/82 cpi transport point.shp")
#cpi <- st_read("Shape file CPI/82 cpi general.shp")

Distance based neighbors - minimum neighbor distance threshold: 338 km

Step 1 Prepare spatial objects for subsequent analysis:

plot(cpi)

plotting the first 10 out of 154 attributes; use max.plot = 154 to plot all

class(cpi)

[1] "sf"         "data.frame"

head(cpi)

Simple feature collection with 6 features and 154 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 113 ymin: -8.69 xmax: 116 ymax: -7.99
geographic CRS: WGS 84
  X.x..R           Province   Regency_Ci        Region   X     Y
1   <NA>               Bali Capital City    Java -Bali 115 -8.69
2   <NA> West Nusa Tenggara Capital City Nusa Tenggara 116 -8.60
3   <NA>          East Java      Regency    Java -Bali 114 -8.21
4   <NA>          East Java      Regency    Java -Bali 114 -8.16
5   <NA>               Bali      Regency    Java -Bali 115 -8.12
6   <NA>          East Java         City    Java -Bali 113 -7.99
      reference IDCity   id       city products mo2014m1 mo2014m2
1 Kota Denpasar   5171 5171   Denpasar  General      109      110
2  Kota Mataram   5271 5271    Mataram  General      111      112
3    Banyuwangi   3510 3510 Banyuwangi  General      111      112
4        Jember   3509 3509     Jember  General      111      111
5      Buleleng   5108 5108  Singaraja  General      115      115
6   Kota Malang   3573 3573     Malang  General      111      111
  mo2014m3 mo2014m4 mo2014m5 mo2014m6 mo2014m7 mo2014m8 mo2014m9
1      110      110      110      110      111      111      112
2      111      111      111      111      112      113      113
3      112      112      112      113      113      113      113
4      111      111      111      111      112      112      112
5      115      115      117      116      117      118      119
6      112      112      112      112      113      114      114
  mo2014m10 mo2014m11 mo2014m12 mo2015m1 mo2015m2 mo2015m3 mo2015m4
1       112       114       116      116      116      116      117
2       114       115       117      118      117      118      118
3       113       115       118      118      117      117      117
4       112       114       118      117      117      117      117
5       120       122       125      125      125      126      126
6       114       116       119      119      119      119      120
  mo2015m5 mo2015m6 mo2015m7 mo2015m8 mo2015m9 mo2015m10 mo2015m11
1      117      117      119      119      119       118       118
2      118      118      119      119      120       120       120
3      118      118      119      119      119       119       119
4      117      118      119      119      120       119       120
5      127      126      128      128      128       127       127
6      120      121      121      122      122       122       122
  mo2015m12 mo2016m1 mo2016m2 mo2016m3 mo2016m4 mo2016m5 mo2016m6
1       120      120      120      120      120      120      121
2       121      123      122      122      122      122      123
3       120      121      121      121      120      121      121
4       120      121      121      121      120      121      121
5       129      131      130      131      131      131      131
6       123      124      124      124      123      123      124
  mo2016m7 mo2016m8 mo2016m9 mo2016m10 mo2016m11 mo2016m12 mo2017m1
1      121      122      122       122       122       123      125
2      124      123      123       123       123       124      126
3      122      122      122       122       122       122      123
4      121      121      121       121       121       123      124
5      132      134      134       133       134       135      138
6      125      125      125       125       126       126      128
  mo2017m2 mo2017m3 mo2017m4 mo2017m5 mo2017m6 mo2017m7 mo2017m8
1      125      125      125      126      126      126      126
2      127      126      126      126      127      128      127
3      124      123      124      124      125      125      125
4      125      124      125      125      126      126      126
5      139      138      137      137      136      137      137
6      128      128      129      130      130      131      130
  mo2017m9 mo2017m10 mo2017m11 mo2017m12 mo2018m1 mo2018m2 mo2018m3
1      126       126       126       127      128      129      129
2      127       128       128       129      129      130      130
3      125       125       126       126      127      127      128
4      126       126       126       127      128      128      128
5      136       136       138       140      141      141      142
6      130       130       130       131      132      132      132
  mo2018m4 mo2018m5 mo2018m6 mo2018m7 mo2018m8 mo2018m9 mo2018m10
1      129      129      130      131      131      130       130
2      130      130      131      132      132      131       132
3      128      128      128      129      128      128       128
4      128      129      130      129      129      129       130
5      141      141      141      142      142      141       141
6      133      133      133      134      134      133       134
  mo2018m11 mo2018m12 mo2019m1 mo2019m2 mo2019m3 mo2019m4 mo2019m5
1       130       132      132      132      132      132      133
2       132       133      133      133      133      133      134
3       128       129      129      129      130      130      130
4       130       131      131      131      131      131      132
5       141       142      143      143      143      144      144
6       134       135      136      135      136      136      137
  mo2019m6 mo2019m7 mo2019m8 mo2019m9 mo2019m10 mo2019m11 mo2019m12
1      133      134      134      133       134       134       135
2      135      135      135      134       135       135       135
3      131      131      131      131       131       132       132
4      132      132      132      132       132       133       133
5      144      146      146      145       145       145       146
6      136      137      137      137       137       137       138
  ln2014m1 ln2014m2 ln2014m3 ln2014m4 ln2014m5 ln2014m6 ln2014m7
1     4.69     4.70     4.70     4.70     4.70     4.70     4.71
2     4.71     4.71     4.71     4.71     4.71     4.71     4.72
3     4.71     4.72     4.72     4.72     4.72     4.72     4.73
4     4.71     4.71     4.71     4.71     4.71     4.71     4.72
5     4.74     4.75     4.75     4.75     4.76     4.76     4.76
6     4.71     4.71     4.72     4.72     4.72     4.72     4.73
  ln2014m8 ln2014m9 ln2014m10 ln2014m11 ln2014m12 ln2015m1 ln2015m2
1     4.71     4.72      4.72      4.74      4.76     4.76     4.76
2     4.73     4.73      4.73      4.74      4.77     4.77     4.77
3     4.72     4.73      4.73      4.74      4.77     4.77     4.76
4     4.72     4.72      4.72      4.74      4.77     4.76     4.76
5     4.77     4.78      4.78      4.80      4.83     4.83     4.83
6     4.73     4.73      4.74      4.75      4.78     4.78     4.78
  ln2015m3 ln2015m4 ln2015m5 ln2015m6 ln2015m7 ln2015m8 ln2015m9
1     4.76     4.76     4.76     4.77     4.77     4.78     4.78
2     4.77     4.77     4.77     4.77     4.78     4.78     4.79
3     4.76     4.76     4.77     4.77     4.78     4.78     4.78
4     4.76     4.76     4.77     4.77     4.78     4.78     4.78
5     4.83     4.84     4.84     4.84     4.85     4.85     4.85
6     4.78     4.78     4.79     4.79     4.80     4.80     4.80
  ln2015m10 ln2015m11 ln2015m12 ln2016m1 ln2016m2 ln2016m3 ln2016m4
1      4.77      4.77      4.78     4.79     4.79     4.79     4.79
2      4.79      4.79      4.80     4.81     4.81     4.81     4.80
3      4.78      4.78      4.79     4.80     4.80     4.80     4.79
4      4.78      4.79      4.79     4.79     4.80     4.80     4.79
5      4.84      4.85      4.86     4.87     4.87     4.88     4.88
6      4.80      4.80      4.81     4.82     4.82     4.82     4.81
  ln2016m5 ln2016m6 ln2016m7 ln2016m8 ln2016m9 ln2016m10 ln2016m11
1     4.79     4.79     4.80     4.80     4.81      4.80      4.81
2     4.80     4.81     4.82     4.82     4.81      4.81      4.82
3     4.79     4.80     4.80     4.80     4.80      4.80      4.80
4     4.79     4.80     4.80     4.80     4.80      4.80      4.80
5     4.88     4.88     4.89     4.89     4.90      4.89      4.90
6     4.82     4.82     4.83     4.83     4.83      4.83      4.83
  ln2016m12 ln2017m1 ln2017m2 ln2017m3 ln2017m4 ln2017m5 ln2017m6
1      4.81     4.83     4.83     4.83     4.83     4.83     4.83
2      4.82     4.84     4.84     4.84     4.83     4.84     4.84
3      4.81     4.81     4.82     4.82     4.82     4.82     4.83
4      4.81     4.82     4.83     4.82     4.83     4.83     4.83
5      4.91     4.92     4.93     4.93     4.92     4.92     4.92
6      4.84     4.85     4.86     4.85     4.86     4.87     4.87
  ln2017m7 ln2017m8 ln2017m9 ln2017m10 ln2017m11 ln2017m12 ln2018m1
1     4.83     4.84     4.83      4.83      4.83      4.85     4.85
2     4.85     4.85     4.85      4.85      4.85      4.86     4.86
3     4.83     4.83     4.83      4.83      4.83      4.84     4.85
4     4.84     4.83     4.83      4.83      4.84      4.84     4.85
5     4.92     4.92     4.91      4.91      4.93      4.94     4.95
6     4.87     4.87     4.87      4.87      4.87      4.88     4.88
  ln2018m2 ln2018m3 ln2018m4 ln2018m5 ln2018m6 ln2018m7 ln2018m8
1     4.86     4.86     4.86     4.86     4.87     4.87     4.87
2     4.87     4.87     4.87     4.87     4.87     4.88     4.88
3     4.85     4.85     4.85     4.85     4.86     4.86     4.86
4     4.85     4.85     4.85     4.86     4.86     4.86     4.86
5     4.95     4.95     4.95     4.95     4.95     4.95     4.96
6     4.88     4.89     4.89     4.89     4.89     4.89     4.90
  ln2018m9 ln2018m10 ln2018m11 ln2018m12 ln2019m1 ln2019m2 ln2019m3
1     4.87      4.87      4.87      4.88     4.89     4.88     4.88
2     4.88      4.88      4.88      4.89     4.89     4.89     4.89
3     4.85      4.85      4.85      4.86     4.86     4.86     4.86
4     4.86      4.86      4.87      4.87     4.87     4.87     4.87
5     4.95      4.95      4.95      4.96     4.96     4.96     4.96
6     4.89      4.89      4.90      4.91     4.91     4.91     4.91
  ln2019m4 ln2019m5 ln2019m6 ln2019m7 ln2019m8 ln2019m9 ln2019m10
1     4.89     4.89     4.89     4.89     4.90     4.89      4.89
2     4.89     4.90     4.91     4.91     4.90     4.90      4.90
3     4.87     4.87     4.87     4.88     4.88     4.88      4.88
4     4.88     4.88     4.88     4.88     4.89     4.88      4.88
5     4.97     4.97     4.97     4.98     4.99     4.98      4.98
6     4.91     4.92     4.92     4.92     4.92     4.92      4.92
  ln2019m11          geometry
1      4.89 POINT (115 -8.69)
2      4.90  POINT (116 -8.6)
3      4.88 POINT (114 -8.21)
4      4.89 POINT (114 -8.16)
5      4.98 POINT (115 -8.12)
6      4.92 POINT (113 -7.99)

tail(cpi)

Simple feature collection with 6 features and 154 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 120 ymin: -0.919 xmax: 134 ymax: 1.56
geographic CRS: WGS 84
   X.x..R         Province   Regency_Ci   Region   X      Y
77   <NA>      Papua Barat Capital City    Papua 131 -0.919
78   <NA>      Papua Barat      Regency    Papua 134 -0.860
79   <NA> Central Sulawesi Capital City Sulawesi 120 -0.808
80   <NA>        Gorontalo Capital City Sulawesi 123  0.547
81   <NA>     North Maluku Capital City   Maluku 127  0.811
82   <NA>   North Sulawesi Capital City Sulawesi 125  1.563
              reference IDCity   id      city products mo2014m1
77          Kota Sorong   9171 9171    Sorong  General      108
78            Manokwari   9105 9105 Manokwari  General      106
79            Kota Palu   7271 7271      Palu  General      112
80 Kota Gorontalo, Kota   7571 7571 Gorontalo  General      109
81   Kota Ternate, Kota   8271 8271   Ternate  General      112
82          Kota Manado   7171 7171    Manado  General      109
   mo2014m2 mo2014m3 mo2014m4 mo2014m5 mo2014m6 mo2014m7 mo2014m8
77      109      109      110      110      110      112      114
78      107      106      106      107      107      108      110
79      111      111      112      113      114      115      116
80      108      108      109      109      109      110      110
81      112      112      113      113      114      117      116
82      109      109      110      110      110      111      111
   mo2014m9 mo2014m10 mo2014m11 mo2014m12 mo2015m1 mo2015m2 mo2015m3
77      115       114       114       116      116      117      117
78      110       111       111       113      112      112      113
79      115       117       117       120      120      118      117
80      110       110       111       115      114      113      114
81      117       118       119       122      122      121      121
82      111       112       114       119      118      118      118
   mo2015m4 mo2015m5 mo2015m6 mo2015m7 mo2015m8 mo2015m9 mo2015m10
77      117      117      120      122      123      123       123
78      113      113      114      115      113      114       113
79      118      120      120      122      121      121       122
80      114      115      116      117      118      118       118
81      122      123      124      125      127      125       126
82      118      119      120      121      121      121       123
   mo2015m11 mo2015m12 mo2016m1 mo2016m2 mo2016m3 mo2016m4 mo2016m5
77       122       123      125      125      125      124      123
78       113       116      116      116      116      116      117
79       123       125      125      124      124      124      125
80       118       120      120      120      120      120      120
81       126       128      128      127      128      128      128
82       123       125      125      124      124      123      123
   mo2016m6 mo2016m7 mo2016m8 mo2016m9 mo2016m10 mo2016m11 mo2016m12
77      124      126      127      127       126       126       127
78      119      120      122      121       120       121       122
79      126      126      126      126       125       126       127
80      122      122      121      121       120       121       122
81      128      130      130      130       130       130       130
82      124      125      125      124       124       128       126
   mo2017m1 mo2017m2 mo2017m3 mo2017m4 mo2017m5 mo2017m6 mo2017m7
77      128      128      129      128      128      129      130
78      122      122      122      121      122      124      125
79      129      129      129      130      131      132      132
80      123      124      124      124      124      126      127
81      131      131      131      131      131      133      135
82      127      128      129      129      127      129      130
   mo2017m8 mo2017m9 mo2017m10 mo2017m11 mo2017m12 mo2018m1 mo2018m2
77      129      129       129       128       129      129      130
78      123      125       124       124       125      126      124
79      132      132       130       130       133      134      133
80      126      126       126       126       127      128      127
81      133      132       133       131       133      134      134
82      130      128       128       128       129      129      130
   mo2018m3 mo2018m4 mo2018m5 mo2018m6 mo2018m7 mo2018m8 mo2018m9
77      130      131      132      134      136      136      135
78      125      125      126      127      128      128      128
79      133      134      134      137      137      137      135
80      127      127      128      129      129      129      129
81      135      136      136      139      137      137      137
82      130      132      132      133      132      131      130
   mo2018m10 mo2018m11 mo2018m12 mo2019m1 mo2019m2 mo2019m3 mo2019m4
77       135       135       135      135      134      133      134
78       129       130       132      133      133      133      133
79       138       140       141      141      141      140      141
80       129       129       130      130      129      129      130
81       137       137       138      139      139      139      139
82       130       133       134      135      134      133      132
   mo2019m5 mo2019m6 mo2019m7 mo2019m8 mo2019m9 mo2019m10 mo2019m11
77      135      136      136      137      137       137       135
78      136      136      135      136      136       137       137
79      143      144      143      144      143       143       143
80      132      132      132      133      133       133       133
81      140      141      141      141      140       140       141
82      135      140      138      136      135       136       141
   mo2019m12 ln2014m1 ln2014m2 ln2014m3 ln2014m4 ln2014m5 ln2014m6
77       136     4.69     4.69     4.69     4.70     4.70     4.70
78       138     4.67     4.67     4.67     4.67     4.67     4.68
79       144     4.71     4.71     4.71     4.72     4.72     4.73
80       134     4.69     4.68     4.68     4.69     4.69     4.69
81       141     4.72     4.71     4.72     4.73     4.73     4.74
82       138     4.69     4.69     4.69     4.70     4.70     4.70
   ln2014m7 ln2014m8 ln2014m9 ln2014m10 ln2014m11 ln2014m12 ln2015m1
77     4.72     4.74     4.75      4.74      4.74      4.75     4.76
78     4.69     4.70     4.70      4.71      4.71      4.72     4.72
79     4.75     4.75     4.75      4.76      4.76      4.79     4.79
80     4.70     4.70     4.70      4.70      4.71      4.75     4.73
81     4.76     4.75     4.76      4.77      4.78      4.81     4.80
82     4.71     4.71     4.71      4.72      4.74      4.78     4.77
   ln2015m2 ln2015m3 ln2015m4 ln2015m5 ln2015m6 ln2015m7 ln2015m8
77     4.76     4.76     4.76     4.77     4.78     4.80     4.81
78     4.72     4.73     4.72     4.72     4.74     4.75     4.73
79     4.77     4.77     4.77     4.79     4.79     4.80     4.80
80     4.73     4.74     4.74     4.75     4.75     4.76     4.77
81     4.79     4.80     4.80     4.81     4.82     4.83     4.84
82     4.77     4.77     4.77     4.78     4.79     4.80     4.79
   ln2015m9 ln2015m10 ln2015m11 ln2015m12 ln2016m1 ln2016m2 ln2016m3
77     4.81      4.81      4.81      4.81     4.82     4.83     4.82
78     4.73      4.73      4.73      4.75     4.75     4.75     4.75
79     4.80      4.81      4.81      4.83     4.83     4.82     4.82
80     4.77      4.77      4.77      4.79     4.78     4.79     4.79
81     4.83      4.84      4.84      4.85     4.86     4.85     4.85
82     4.80      4.81      4.81      4.83     4.83     4.82     4.82
   ln2016m4 ln2016m5 ln2016m6 ln2016m7 ln2016m8 ln2016m9 ln2016m10
77     4.82     4.81     4.82     4.83     4.85     4.85      4.84
78     4.75     4.76     4.78     4.79     4.80     4.79      4.79
79     4.82     4.83     4.83     4.84     4.83     4.84      4.83
80     4.79     4.79     4.80     4.80     4.80     4.80      4.79
81     4.85     4.85     4.86     4.87     4.86     4.87      4.86
82     4.81     4.81     4.82     4.83     4.83     4.82      4.82
   ln2016m11 ln2016m12 ln2017m1 ln2017m2 ln2017m3 ln2017m4 ln2017m5
77      4.84      4.84     4.85     4.85     4.86     4.85     4.85
78      4.80      4.81     4.81     4.80     4.80     4.80     4.81
79      4.83      4.84     4.86     4.86     4.86     4.87     4.88
80      4.80      4.80     4.81     4.82     4.82     4.82     4.82
81      4.87      4.87     4.88     4.88     4.87     4.88     4.88
82      4.85      4.83     4.84     4.86     4.86     4.86     4.85
   ln2017m6 ln2017m7 ln2017m8 ln2017m9 ln2017m10 ln2017m11 ln2017m12
77     4.86     4.86     4.86     4.86      4.86      4.85      4.86
78     4.82     4.83     4.81     4.83      4.82      4.82      4.82
79     4.88     4.88     4.88     4.88      4.87      4.87      4.89
80     4.84     4.85     4.84     4.84      4.84      4.84      4.84
81     4.89     4.90     4.89     4.88      4.89      4.88      4.89
82     4.86     4.87     4.86     4.85      4.85      4.85      4.86
   ln2018m1 ln2018m2 ln2018m3 ln2018m4 ln2018m5 ln2018m6 ln2018m7
77     4.86     4.87     4.87     4.88     4.88     4.90     4.91
78     4.83     4.82     4.82     4.83     4.84     4.85     4.85
79     4.89     4.89     4.89     4.90     4.90     4.92     4.92
80     4.85     4.84     4.85     4.85     4.85     4.86     4.86
81     4.90     4.90     4.91     4.91     4.92     4.93     4.92
82     4.86     4.87     4.87     4.88     4.89     4.89     4.89
   ln2018m8 ln2018m9 ln2018m10 ln2018m11 ln2018m12 ln2019m1 ln2019m2
77     4.92     4.90      4.90      4.91      4.90     4.91     4.90
78     4.85     4.85      4.86      4.87      4.88     4.89     4.89
79     4.92     4.91      4.93      4.94      4.95     4.95     4.95
80     4.86     4.86      4.86      4.86      4.87     4.87     4.86
81     4.92     4.92      4.92      4.92      4.93     4.94     4.93
82     4.88     4.87      4.87      4.89      4.90     4.91     4.90
   ln2019m3 ln2019m4 ln2019m5 ln2019m6 ln2019m7 ln2019m8 ln2019m9
77     4.89     4.89     4.91     4.91     4.92     4.92     4.92
78     4.89     4.89     4.91     4.91     4.91     4.92     4.91
79     4.94     4.95     4.96     4.97     4.96     4.97     4.96
80     4.86     4.87     4.88     4.89     4.89     4.89     4.89
81     4.93     4.94     4.94     4.95     4.95     4.95     4.94
82     4.89     4.88     4.91     4.94     4.93     4.91     4.90
   ln2019m10 ln2019m11           geometry
77      4.92      4.91 POINT (131 -0.919)
78      4.92      4.92  POINT (134 -0.86)
79      4.96      4.96 POINT (120 -0.808)
80      4.89      4.89  POINT (123 0.547)
81      4.94      4.95  POINT (127 0.811)
82      4.92      4.95   POINT (125 1.56)

#minimum distance is 338 for the dnearneigh fuction

#convert cpi spatial object to data.frame

cpi.df <- as.data.frame(cpi)

nb <- dnearneigh(cpi, d1=0, d2=672+ 10*jj, longlat=T, row.names = IDs) #create weight matrix and name W.matrix

# (a) coordinates and IDs
#coords <- cpi[,c("X", "Y")]
#coords <- coordinates(coords)
IDs <- cpi$IDCity
# (b) identify neighbors given the distance threshold
nb672km <- dnearneigh(cpi, d1=0, d2=672, row.names = IDs)
summary(nb672km)

Neighbour list object:
Number of regions: 82 
Number of nonzero links: 1588 
Percentage nonzero weights: 23.6 
Average number of links: 19.4 
Link number distribution:

 1  2  3  4  5  6  7  9 10 12 13 14 15 16 17 18 19 20 23 24 25 
 2  1  1  5  1  1  2  2  2  1  3  2  3  7  1  5  4  2  3  2  3 
26 27 28 29 30 31 32 33 34 
 3  5  3  8  2  3  2  1  2 
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links

# (c) calculate the spatial weights matrix
W.matrixdist <- nb2listw(nb672km)
summary(W.matrixdist)

Characteristics of weights list object:
Neighbour list object:
Number of regions: 82 
Number of nonzero links: 1588 
Percentage nonzero weights: 23.6 
Average number of links: 19.4 
Link number distribution:

 1  2  3  4  5  6  7  9 10 12 13 14 15 16 17 18 19 20 23 24 25 
 2  1  1  5  1  1  2  2  2  1  3  2  3  7  1  5  4  2  3  2  3 
26 27 28 29 30 31 32 33 34 
 3  5  3  8  2  3  2  1  2 
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links

Weights style: W 
Weights constants summary:

#nb2listw(nb338km) #, zero.policy = TRUE)
#W.matrix <- nb2listw(nb338km,zero.policy = TRUE)
#summary(W.matrix)

#weight matrix k nearest

# (a) coordinates and IDs
#coords <- cpi[,c("X", "Y")]
#coords <- coordinates(coords)
IDs <- cpi$IDCity
# (b) identify neighbors given the distance threshold
nearneigh4 <- knearneigh(cpi, k=4)
summary(nearneigh4)

          Length Class  Mode   
nn        328    -none- numeric
np          1    -none- numeric
k           1    -none- numeric
dimension   1    -none- numeric
x         164    -none- numeric

# (c) calculate the spatial weights matrix
W.matrixnearest <- nb2listw(knn2nb(nearneigh4))
summary(W.matrixnearest)

Characteristics of weights list object:
Neighbour list object:
Number of regions: 82 
Number of nonzero links: 328 
Percentage nonzero weights: 4.88 
Average number of links: 4 
Non-symmetric neighbours list
Link number distribution:

 4 
82 
82 least connected regions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 with 4 links
82 most connected regions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 with 4 links

Weights style: W 
Weights constants summary:

#checking the average of Moran’s I

cpi <- cpi %>% 
  mutate(avg = (mo2016m10+mo2016m11+mo2016m12+mo2017m1+mo2017m2+mo2017m3)/6)

Step 2 Spatial filtering

Getis filtering can be applied for

positively autocorrelated data with natural origin:

#compute Moran’s I for each observation

#moran.test(cpi$X2014m1,W.matrix)
#moran.test(cpi$X2015m1,W.matrix)
#moran.test(cpi$X2016m1,W.matrix)
#moran.test(cpi$X2017m1,W.matrix)
#moran.test(cpi$X2018m1,W.matrix)
#moran.test(cpi$X2019m1,W.matrix)

moran.test(cpi$mo2016m10,W.matrix)

    Moran I test under randomisation

data:  cpi$mo2016m10  
weights: W.matrix    

Moran I statistic standard deviate = 0.7, p-value = 0.2
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
          0.01932          -0.01235           0.00203 

moran.test(cpi$mo2016m12,W.matrix)

    Moran I test under randomisation

data:  cpi$mo2016m12  
weights: W.matrix    

Moran I statistic standard deviate = 1, p-value = 0.1
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
          0.04559          -0.01235           0.00202 

moran.test(cpi$mo2017m2,W.matrix)

    Moran I test under randomisation

data:  cpi$mo2017m2  
weights: W.matrix    

Moran I statistic standard deviate = -0.5, p-value = 0.7
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
         -0.03265          -0.01235           0.00202 

moran.test(cpi$avg,W.matrix)

    Moran I test under randomisation

data:  cpi$avg  
weights: W.matrix    

Moran I statistic standard deviate = 0.5, p-value = 0.3
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
          0.01163          -0.01235           0.00202 

names=c()
names <- colnames(cpi.df)
Moran.s_I=data.frame(var="a",moran=0, p.value=0)
for (i in 12:83) {
    MOR<- moran.test(cpi.df[,i],W.matrixnearest)
  lab= (round(MOR$estimate[1],3))
  pvalue = MOR$p.value
a=data.frame(var=names[i],moran=lab, p.value=pvalue)
Moran.s_I<-rbind(Moran.s_I,a)
}
Moran.s_I<-Moran.s_I[-1,]
tail(Moran.s_I)

Moran.s_I

filter(Moran.s_I,p.value < 0.1)

Spatial autocorrelation

General = 28 out of 72 (k=4)

Foodstuff = 52

Procfood = 8

Housing = 35

Health = 63

Clothing = 0

Education = 0

Transport = 6

cpi.df 

cpi[,11]

Simple feature collection with 82 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 95.3 ymin: -10.2 xmax: 141 ymax: 5.56
CRS:            4326
First 10 features:
   mo2014m1          geometry
1       116 POINT (115 -8.69)
2       120  POINT (116 -8.6)
3       127 POINT (114 -8.21)
4       117 POINT (114 -8.16)
5       118 POINT (115 -8.12)
6       118 POINT (113 -7.99)
7       116 POINT (112 -7.83)
8       120 POINT (110 -7.79)
9       117 POINT (113 -7.75)
10      120  POINT (109 -7.7)

Testing the Getis G

# (b) identify neighbors given the distance threshold
nb672km <- dnearneigh(cpi, d1=0, d2=672, row.names = IDs)
summary(nb672km)

Neighbour list object:
Number of regions: 82 
Number of nonzero links: 1588 
Percentage nonzero weights: 23.6 
Average number of links: 19.4 
Link number distribution:

 1  2  3  4  5  6  7  9 10 12 13 14 15 16 17 18 19 20 23 24 25 
 2  1  1  5  1  1  2  2  2  1  3  2  3  7  1  5  4  2  3  2  3 
26 27 28 29 30 31 32 33 34 
 3  5  3  8  2  3  2  1  2 
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links

# (c) calculate the spatial weights matrix
W.matrix <- nb2listw(nb672km)
summary(W.matrix)

Characteristics of weights list object:
Neighbour list object:
Number of regions: 82 
Number of nonzero links: 1588 
Percentage nonzero weights: 23.6 
Average number of links: 19.4 
Link number distribution:

 1  2  3  4  5  6  7  9 10 12 13 14 15 16 17 18 19 20 23 24 25 
 2  1  1  5  1  1  2  2  2  1  3  2  3  7  1  5  4  2  3  2  3 
26 27 28 29 30 31 32 33 34 
 3  5  3  8  2  3  2  1  2 
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links

Weights style: W 
Weights constants summary:

?localG
localG(cpi.df$mo2014m3, W.matrix)

 [1]  1.64878  1.38398  1.15990  1.85607  1.38404  2.23471
 [7]  2.74276  4.20782  1.88619  4.18377  2.90712  4.10851
[13]  4.11783  2.96376  2.32004  1.54396  3.87443  3.04638
[19]  2.84404  3.55592  3.73188  3.44055  2.83751  2.74521
[25]  2.74096  2.59241  2.64536  2.66724  2.59283  3.15398
[31]  2.79039  2.11122  0.33211  2.41516  2.63926  3.27713
[37]  1.20955 -1.13287  1.12867 -2.04853 -1.28549 -0.92493
[43] -2.56634 -3.12814 -1.45965 -3.13347 -2.42568 -2.84249
[49] -2.48415 -3.21206 -2.68763  0.34423  0.81283  2.67875
[55] -1.31488  2.28088  1.37271 -1.58224 -0.27780 -2.16889
[61]  0.63526 -0.22797 -0.57143 -0.00653  0.28373 -0.05982
[67]  0.06157 -0.27498 -0.73210 -0.32900 -0.30385 -0.64036
[73] -0.35039 -1.28095 -1.35349 -1.37972 -1.34818 -0.83486
[79] -2.45196 -2.64425 -1.16874 -3.28812
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cpi.df$mo2014m3, listw = W.matrix)
attr(,"class")
[1] "localG"

getis<- localG(cpi.df$mo2014m3, nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
getis

 [1]  1.64878  1.38398  1.15990  1.85607  1.38404  2.23471
 [7]  2.74276  4.20782  1.88619  4.18377  2.90712  4.10851
[13]  4.11783  2.96376  2.32004  1.54396  3.87443  3.04638
[19]  2.84404  3.55592  3.73188  3.44055  2.83751  2.74521
[25]  2.74096  2.59241  2.64536  2.66724  2.59283  3.15398
[31]  2.79039  2.11122  0.33211  2.41516  2.63926  3.27713
[37]  1.20955 -1.13287  1.12867 -2.04853 -1.28549 -0.92493
[43] -2.56634 -3.12814 -1.45965 -3.13347 -2.42568 -2.84249
[49] -2.48415 -3.21206 -2.68763  0.34423  0.81283  2.67875
[55] -1.31488  2.28088  1.37271 -1.58224 -0.27780 -2.16889
[61]  0.63526 -0.22797 -0.57143 -0.00653  0.28373 -0.05982
[67]  0.06157 -0.27498 -0.73210 -0.32900 -0.30385 -0.64036
[73] -0.35039 -1.28095 -1.35349 -1.37972 -1.34818 -0.83486
[79] -2.45196 -2.64425 -1.16874 -3.28812
attr(,"internals")
           G     EG          VG
 [1,] 0.2261 0.2222 0.000005453
 [2,] 0.2503 0.2469 0.000005852
 [3,] 0.2868 0.2840 0.000005949
 [4,] 0.2887 0.2840 0.000006415
 [5,] 0.2875 0.2840 0.000006406
 [6,] 0.3392 0.3333 0.000007000
 [7,] 0.3654 0.3580 0.000007245
 [8,] 0.4066 0.3951 0.000007499
 [9,] 0.3011 0.2963 0.000006579
[10,] 0.3693 0.3580 0.000007242
[11,] 0.3907 0.3827 0.000007454
[12,] 0.4187 0.4074 0.000007611
[13,] 0.3815 0.3704 0.000007336
[14,] 0.3660 0.3580 0.000007247
[15,] 0.3767 0.3704 0.000007347
[16,] 0.3126 0.3086 0.000006730
[17,] 0.4305 0.4198 0.000007681
[18,] 0.3538 0.3457 0.000007129
[19,] 0.3533 0.3457 0.000007128
[20,] 0.3924 0.3827 0.000007441
[21,] 0.4297 0.4198 0.000007071
[22,] 0.3921 0.3827 0.000007435
[23,] 0.3657 0.3580 0.000007252
[24,] 0.3654 0.3580 0.000007191
[25,] 0.3654 0.3580 0.000007185
[26,] 0.3402 0.3333 0.000006951
[27,] 0.3527 0.3457 0.000007047
[28,] 0.3279 0.3210 0.000006811
[29,] 0.3402 0.3333 0.000007012
[30,] 0.3417 0.3333 0.000007006
[31,] 0.3400 0.3333 0.000005643
[32,] 0.3017 0.2963 0.000006553
[33,] 0.3219 0.3210 0.000006848
[34,] 0.3149 0.3086 0.000006674
[35,] 0.3651 0.3580 0.000007125
[36,] 0.4040 0.3951 0.000007473
[37,] 0.3118 0.3086 0.000006732
[38,] 0.0480 0.0494 0.000001459
[39,] 0.0130 0.0123 0.000000384
[40,] 0.1076 0.1111 0.000002996
[41,] 0.2192 0.2222 0.000005450
[42,] 0.0361 0.0370 0.000001124
[43,] 0.2162 0.2222 0.000005454
[44,] 0.1540 0.1605 0.000004245
[45,] 0.2311 0.2346 0.000005600
[46,] 0.1905 0.1975 0.000004994
[47,] 0.1557 0.1605 0.000003948
[48,] 0.1912 0.1975 0.000004937
[49,] 0.0584 0.0617 0.000001819
[50,] 0.1782 0.1852 0.000004757
[51,] 0.1916 0.1975 0.000004823
[52,] 0.0126 0.0123 0.000000379
[53,] 0.2365 0.2346 0.000005588
[54,] 0.3652 0.3580 0.000007077
[55,] 0.1946 0.1975 0.000004996
[56,] 0.3270 0.3210 0.000006850
[57,] 0.2378 0.2346 0.000005625
[58,] 0.1572 0.1605 0.000004240
[59,] 0.1969 0.1975 0.000004981
[60,] 0.1682 0.1728 0.000004511
[61,] 0.2484 0.2469 0.000005539
[62,] 0.1970 0.1975 0.000004994
[63,] 0.1101 0.1111 0.000003116
[64,] 0.1975 0.1975 0.000004976
[65,] 0.1240 0.1235 0.000003413
[66,] 0.2221 0.2222 0.000005329
[67,] 0.2347 0.2346 0.000005662
[68,] 0.2092 0.2099 0.000005206
[69,] 0.2205 0.2222 0.000005454
[70,] 0.1845 0.1852 0.000004758
[71,] 0.1722 0.1728 0.000004508
[72,] 0.0486 0.0494 0.000001476
[73,] 0.1474 0.1481 0.000003971
[74,] 0.0844 0.0864 0.000002439
[75,] 0.0843 0.0864 0.000002455
[76,] 0.0721 0.0741 0.000002123
[77,] 0.0478 0.0494 0.000001465
[78,] 0.0240 0.0247 0.000000741
[79,] 0.1798 0.1852 0.000004735
[80,] 0.1187 0.1235 0.000003268
[81,] 0.0480 0.0494 0.000001427
[82,] 0.0454 0.0494 0.000001481
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cpi.df$mo2014m3, listw = nb2listw(nb672km, style = "B"), 
    return_internals = T, GeoDa = T)
attr(,"class")
[1] "localG"

For spatial filtering, we need

- return_internals = TRUE

.. G(i) and E(G) are used in filtering

- GeoDa=FALSE = T

.. exclude “self-neighborhood” - use G(i) not G*(i)

- filtering is usually calculated based on connectivity matrix (C, not W)

.. see syntax below

a<- localG(dist$RGDPPC2017, nb2listw(nb338km, style="B"), return_internals = T, GeoDa = T)
#a
a[1]
#attributes(a)
b<-attributes(a)$internals
#b

Spatial filter for GDPpc years 00,01,01, 05, 09, 13, 15, 16 and 17 performed at distance 338 km

mo2014m1.g <- localG(cpi.df$mo2014m1, nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
Getis.m2014m1 <- as.data.frame(attr(mo2014m1.g, which = "internals")) # retrieve "internals"
cpi.df$g.2014m1 <- cpi.df$mo2014m1 * (Getis.m2014m1$EG/Getis.m2014m1$G ) # "multiplicative filter"

filtered <- list()

a=c()
a<- cpi.df[,7]

a

 [1] 5171 5271 3510 3509 5108 3573 3571 3471 3574 3301 3577
[12] 3372 3302 3278 3578 3529 3374 3272 3273 3376 3319 3274
[23] 3271 3276 3275 3671 3173 3673 3672 1871 1872 1771 6371
[34] 1674 1671 1902 6202 5371 9401 5310 5272 8172 7302 7472
[45] 7371 7311 7471 7372 8171 7373 7604 9471 6271 1971 6309
[56] 1571 1509 6471 1371 6472 1403 1375 6171 1471 6172 2172
[67] 2171 1277 1473 1271 1273 6571 1275 1107 1174 1171 9171
[78] 9105 7271 7571 8271 7171

for (i in 12:82) {
  g <- localG(cpi.df[,i], nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
Getis <- as.data.frame(attr(g, which = "internals")) 
filtered[[i]] <- cpi.df[,i] * (Getis$EG/Getis$G ) # "multiplicative filter"
  a<- as.data.frame(cbind(a,filtered[[i]])) 
  colnames(a)[ncol(a)] <- paste("f",colnames(cpi.df)[i],sep = "_")
}

a

combined_cpi_g <- left_join(cpi.df,a, by = c("IDCity" ="a"))

combined_cpi_g

#saving the filtered data in a new file

dist

dist
out<-dist %>% 
  select(1:4, 45,46,47, 50, 54, 58, 60,61,62, 68:85)
out
write.csv(out, "01-raw-data/output_spatial_fitlering/2017filt.csv")

Evaluation of robustness of delta (338 km) when the distances changes

We crearte an “empty” data.frame to collect estimation results at different

max. neighbor distance thresholds

s2.df <- data.frame( V1=0, cpigetis= 0)
s2.df

NA

Start Getis filter

Main calculation: neighbors & W.matrices are calculated for

distance thresholds from 339 km to 839 km (-km iterations).

Model data are stored into the “s2.df” data.frame

.. calculation may take a few moments

#create a data frame of the same type as the one created by the PPA program # please be aware that the code is much longer and complicated than needed but it was better to use what i have done in the past to study the output of the PPA program

W.matrixnearest <- nb2listw(knn2nb(nearneigh4))

cpigetis <- localG(cpi.df$mo2014m3, nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
cpigetis

 [1]  1.64878  1.38398  1.15990  1.85607  1.38404  2.23471
 [7]  2.74276  4.20782  1.88619  4.18377  2.90712  4.10851
[13]  4.11783  2.96376  2.32004  1.54396  3.87443  3.04638
[19]  2.84404  3.55592  3.73188  3.44055  2.83751  2.74521
[25]  2.74096  2.59241  2.64536  2.66724  2.59283  3.15398
[31]  2.79039  2.11122  0.33211  2.41516  2.63926  3.27713
[37]  1.20955 -1.13287  1.12867 -2.04853 -1.28549 -0.92493
[43] -2.56634 -3.12814 -1.45965 -3.13347 -2.42568 -2.84249
[49] -2.48415 -3.21206 -2.68763  0.34423  0.81283  2.67875
[55] -1.31488  2.28088  1.37271 -1.58224 -0.27780 -2.16889
[61]  0.63526 -0.22797 -0.57143 -0.00653  0.28373 -0.05982
[67]  0.06157 -0.27498 -0.73210 -0.32900 -0.30385 -0.64036
[73] -0.35039 -1.28095 -1.35349 -1.37972 -1.34818 -0.83486
[79] -2.45196 -2.64425 -1.16874 -3.28812
attr(,"internals")
           G     EG          VG
 [1,] 0.2261 0.2222 0.000005453
 [2,] 0.2503 0.2469 0.000005852
 [3,] 0.2868 0.2840 0.000005949
 [4,] 0.2887 0.2840 0.000006415
 [5,] 0.2875 0.2840 0.000006406
 [6,] 0.3392 0.3333 0.000007000
 [7,] 0.3654 0.3580 0.000007245
 [8,] 0.4066 0.3951 0.000007499
 [9,] 0.3011 0.2963 0.000006579
[10,] 0.3693 0.3580 0.000007242
[11,] 0.3907 0.3827 0.000007454
[12,] 0.4187 0.4074 0.000007611
[13,] 0.3815 0.3704 0.000007336
[14,] 0.3660 0.3580 0.000007247
[15,] 0.3767 0.3704 0.000007347
[16,] 0.3126 0.3086 0.000006730
[17,] 0.4305 0.4198 0.000007681
[18,] 0.3538 0.3457 0.000007129
[19,] 0.3533 0.3457 0.000007128
[20,] 0.3924 0.3827 0.000007441
[21,] 0.4297 0.4198 0.000007071
[22,] 0.3921 0.3827 0.000007435
[23,] 0.3657 0.3580 0.000007252
[24,] 0.3654 0.3580 0.000007191
[25,] 0.3654 0.3580 0.000007185
[26,] 0.3402 0.3333 0.000006951
[27,] 0.3527 0.3457 0.000007047
[28,] 0.3279 0.3210 0.000006811
[29,] 0.3402 0.3333 0.000007012
[30,] 0.3417 0.3333 0.000007006
[31,] 0.3400 0.3333 0.000005643
[32,] 0.3017 0.2963 0.000006553
[33,] 0.3219 0.3210 0.000006848
[34,] 0.3149 0.3086 0.000006674
[35,] 0.3651 0.3580 0.000007125
[36,] 0.4040 0.3951 0.000007473
[37,] 0.3118 0.3086 0.000006732
[38,] 0.0480 0.0494 0.000001459
[39,] 0.0130 0.0123 0.000000384
[40,] 0.1076 0.1111 0.000002996
[41,] 0.2192 0.2222 0.000005450
[42,] 0.0361 0.0370 0.000001124
[43,] 0.2162 0.2222 0.000005454
[44,] 0.1540 0.1605 0.000004245
[45,] 0.2311 0.2346 0.000005600
[46,] 0.1905 0.1975 0.000004994
[47,] 0.1557 0.1605 0.000003948
[48,] 0.1912 0.1975 0.000004937
[49,] 0.0584 0.0617 0.000001819
[50,] 0.1782 0.1852 0.000004757
[51,] 0.1916 0.1975 0.000004823
[52,] 0.0126 0.0123 0.000000379
[53,] 0.2365 0.2346 0.000005588
[54,] 0.3652 0.3580 0.000007077
[55,] 0.1946 0.1975 0.000004996
[56,] 0.3270 0.3210 0.000006850
[57,] 0.2378 0.2346 0.000005625
[58,] 0.1572 0.1605 0.000004240
[59,] 0.1969 0.1975 0.000004981
[60,] 0.1682 0.1728 0.000004511
[61,] 0.2484 0.2469 0.000005539
[62,] 0.1970 0.1975 0.000004994
[63,] 0.1101 0.1111 0.000003116
[64,] 0.1975 0.1975 0.000004976
[65,] 0.1240 0.1235 0.000003413
[66,] 0.2221 0.2222 0.000005329
[67,] 0.2347 0.2346 0.000005662
[68,] 0.2092 0.2099 0.000005206
[69,] 0.2205 0.2222 0.000005454
[70,] 0.1845 0.1852 0.000004758
[71,] 0.1722 0.1728 0.000004508
[72,] 0.0486 0.0494 0.000001476
[73,] 0.1474 0.1481 0.000003971
[74,] 0.0844 0.0864 0.000002439
[75,] 0.0843 0.0864 0.000002455
[76,] 0.0721 0.0741 0.000002123
[77,] 0.0478 0.0494 0.000001465
[78,] 0.0240 0.0247 0.000000741
[79,] 0.1798 0.1852 0.000004735
[80,] 0.1187 0.1235 0.000003268
[81,] 0.0480 0.0494 0.000001427
[82,] 0.0454 0.0494 0.000001481
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cpi.df$mo2014m3, listw = nb2listw(nb672km, style = "B"), 
    return_internals = T, GeoDa = T)
attr(,"class")
[1] "localG"

s2.df <- data.frame( V1=0, cpigetis= 0)
s2.df

for(jj in 1:30) {
  nb <- dnearneigh(cpi, d1=0, d2=662+ 10*jj, row.names = IDs)
  # Spatial filter for log.GDPPCPCXX or for RGDPPC20XX
  cpigetis <- localG(cpi.df$mo2014m3, nb2listw(nb, style="B"), return_internals = T, GeoDa = T)
  
  #the output is created so that it resembles the output of the PPA program
   aj <- cbind(cpi.df$IDCity, cpigetis ) 
   aj<- as.data.frame(aj)
   aj
   n= c(NA, NA)
   m= c(jj, NA)
  nm<-  data.frame( n,m )
  names(nm)[1] <- "V1"
names(nm)[2] <- "cpigetis"
  aj <- rbind( nm, aj)
  s2.df<- rbind(s2.df, aj)
  
} 
s2.df <- s2.df[-1,]
#
head(s2.df)

tail(s2.df)

#

s2.df

what follows is the same as the R program I created for evaluating delta for the output file of the PPA program

locg<- s2.df%>% 
  select(1,2)

locg

names(locg)[1] <- "point"
names(locg)[2] <- "getis"
locg

in the following code insert the number of points (regions, 514 districs) and the number of different distances used for calculating the local Gi statistic (100 as written int he for loop on line 228 of this code “for(jj in 1:100) {”…)

nump=82
numd=30

rows= (nump+2)*numd
rows

[1] 2520

As you can see this data is very messy, there are two many rows we need two add a new coloumn (r) which refers to the distance used for calulating Gi.

note: just the absolute value of Gi is needed

r<- c(1:rows)

for (n in seq_along(r)) {
  for (j in 1:numd){
if ((nump+2)*(j-1)< n & n<= (nump+2)*j ) {
  r[n]=662+ 10*j
}
} 
}

df <- data.frame(cbind(locg, r))

df<- df %>% 
  select(point,getis,r) %>% 
  filter(!is.na(point)) %>%
  mutate(getis= sqrt(getis*getis)) #just the absolute value of Gi is needed

df

NA

now having the new column r is possible to spread the dataset so that each observations is the getis statictic each distance is a row and each column represents each point

df <- df %>% 
  select(point,getis,r) %>% 
spread (point,getis) 

df

the distance delta is the distance for which the statistic Gi starts to decrease in absolute value

df

dfmat<- as.matrix(df)

dif= matrix(data=1000, nrow = numd-1, ncol = nump+1)

for (i in 1:(numd-10)) {
  for (j in 1:(nump+1)) {
    if ( dfmat[i+1,j]-dfmat[i,j]< -0.01  & dfmat[i+2,j]-dfmat[i,j]< -0.01 & dfmat[i+3,j]-dfmat[i,j] < -0.01 & dfmat[i+4,j]-dfmat[i,j] < -0.01 ) {
    dif[i,j]<- dfmat[i,1]
    }
     }
}
#deleted this line in the if condition  && dfmat[i+2,j]-dfmat[i+1,j] < -0.1
#&& dfmat[i+2,j]-dfmat[i+1,j] <= 0
difframe<-as.data.frame(dif)
 
difframe

NA

the cells with values equal to 1000 do not have any meaning the other cells represent the distances for which the statistic decreases in absolute value.

for each column there may be many distances for which the statitic decreases in value. Just the first of such values is needed.

the critic value of lambda is at the end of the code (the mode of the distances distribution)

mind<- c(1:nump)


for (i in seq_along(mind)) {
 mind[i]= min(dplyr::pull(difframe, i+1)) #Just the first of such values is needed 
}

point<- c(1:nump)
hist2<- data.frame(cbind(point, mind))
#hist2

p<- hist2 %>% 
  filter(mind<1000) %>% 
ggplot()+
geom_histogram(mapping = aes(x = mind), binwidth = 10)

p

library(plotly)
ggplotly(p)

hist2 %>% 
  filter(mind<1000) %>% 
  summarise(mean= mean(mind), median=median(mind), mode=mfv(mind))

NA
NA
NA

#delta is the mode = 338 km

END

sessionInfo()